Transport Bifurcation Induced by Sheared Toroidal Flow in Tokamak Plasmas 
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First-principles numerical simulations are used to describe a transport bifurcation in a differentially rotating 
tokamak plasma. Such a bifurcation is more probable in a region of zero magnetic shear than one of finite 
magnetic shear because in the former case the component of the sheared toroidal flow that is perpendicular to 
the magnetic field has the strongest suppressing effect on the turbulence. In the zero-magnctic-shcar regime, 
there are no growing linear eigenmodcs at any finite value of flow shear. However, subcritical turbulence can 
be sustained, owing to the existence of modes, driven by the ion temperature gradient and the parallel velocity 
gradient, which grow transiently. Nonetheless, in a parameter space containing a wide range of temperature 
gradients and velocity shears, there is a sizeable window where all turbulence is suppressed. Combined 
with the relatively low transport of momentum by coUisional (neoclassical) mechanisms, this produces the 
conditions for a bifurcation from low to high temperature and velocity gradients. A parametric model is 
constructed which accurately describes the combined effect of the temperature gradient and the flow gradient 
over a wide range of their values. Using this parametric model, it is shown that in this reduced-transport 
state, heat is transported almost neoclassically, while momentum transport is dominated by subcritical PVG 
turbulence. It is further shown that for any given input of torque, there is an optimum input of heat which 
maximises the temperature gradient. The parametric model describes both the behaviour of the subcritical 
turbulence (which cannot be modelled by the quasi-linear methods used in current transport codes) and the 
complicated effect of the flow shear on the transport stiffness. It may prove useful for transport modelling of 
tokamaks with sheared flows. 



I. INTRODUCTION 

A major obstacle to the development of a working mag- 
netic confinement fusion device is the transport of heat 
by pressure-gradient driven turbulence^. Experimental 
results indicate that a shear in the equilibrium plasma 
flow can reduce, or even eliminate, this turbulence^"— . 
Several numerical studies (gyrofluid^, quasilinear— , and 
gyrokinetio^"— ) have reproduced this effect and estab- 
lished that it results from the shear in the flow veloc- 
ity component perpendicular to the equilibrium magnetic 
fleld. It has also been noted that a shear in the flow veloc- 
ity component parallel to the field, which leads to a linear 
instabilityi^"— , can have the opposite effect and increase 
the turbulence amplitudes^i^rJi. Some numerical studies 
have concluded that, in the regimes considered, when the 
effect of parallel flow shear is included the turbulence is 
never fully suppressed^ii. This contrasts with experimen- 
tal results where almost neoclassical, i.e., collisionali^, 
levels of heat transport have been observed^i^. 

The ratio between the destabilising parallel and stabil- 
ising perpendicular components of the flow shear is set by 
the angle of the flow with respect to the magnetic field. 
In general, friction between trapped and passing particles 
and magnetic pumping keep the equilibrium flow nearly 
toroidal, with a large component of destabilising paral- 
lel flow shear—"—. Nonetheless, two recent paperaiSili 
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which investigated values of the flow shear several times 
larger than the linear growth rate of the ion temperature 
gradient (ITG) instability (one of the principal drivers of 
turbulence in fusion devices) found that for relatively low 
temperature gradients there is a range of flow shear val- 
ues where a purely toroidal flow can completely quench 
the turbulence. 

The question remains how large shears in the equilib- 
rium toroidal flow can be achieved, given finite available 
sources of angular momentum. In Refs.— the turbu- 
lent flux of angular momentum was calculated and found 
to be large at moderate flow shears; if it is too large, 
increasing the input of angular momentum will not gen- 
erate strong enough shear in the flow. RefJ^ considered 
the standard Cyclone Base Casein, with the normalised 
inverse magnetic gradient scale length s = 0.8 and a 
range of values of the temperature gradient and the ve- 
locity shear. In this paper, which expands the results 
reported inii, we consider a similar regime to that de- 
scribed ini^, with the difference that here the magnetic 
shear s = 0. This choice is motivated by the experi- 
mental studies which have found that internal transport 
barriers (local regions of very steep flow and temperature 
gradients) tend to form in regions of zero or low magnetic 
shear— and by previous numerical results which show 
a lower ion thermal diffusivity at zero magnetic shear—. 
The results reported below show that the range of flow 
velocity gradients and ion temperature gradients where 
the turbulence is suppressed is much larger than in the 
case of flnite magnetic shear and, therefore, that a transi- 
tion in both the flow and temperature gradients from low 
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to high values is more readily achieved (such a transition 
is in principle also possible with s = 0.8, but in a much 
smaller region of parameter space^^). A positive feed- 
back between the increase in the flow gradient and the 
suppression of the turbulence provides the mechanism for 
a jump from low to high flow gradients, with a simultane- 
ous jump in the temperature gradient. We show that this 
transition results in a state where the transport of heat 
is nearly neoclassical, whilst the transport of momentum 
remains largely turbulent. 

In a fusion device, the quantities that can be controlled 
are the input of heat and toroidal angular momentum, 
which in steady state are proportional to their outgo- 
ing fluxes. By varying these quantities, it is possible to 
vary the equilibrium gradients. In contrast, in the local 
gyrokinetic simulations used in this study the input pa- 
rameters are the gradients, and the fluxes are calculated 
from the output. In order to demonstrate the existence of 
a bifurcation in the gradients at fixed values of the input 
fluxes, we use local gyrokinetic simulations to map out 
the dependence of the turbulent fluxes on the tempera- 
ture and flow gradients, and then we add the neoclassical 
fluxes and invert this numerically determined dependence 
to find the gradients as functions of the total fluxes. We 
first do this by straightforward interpolation in the pa- 
rameter space, then propose a simple parameterisation of 
the fluxes that fits the data and can be used to predict 
in which parameter regimes bifurcations can occur. 

The rest of this paper is organised as follows. In Sec- 
tion ini we present the gyrokinetic system of equations 
used in our simulations and describe the numerical setup. 
In Section [nil we report a numerical scan in two param- 
eters: the flow shear and the ion temperature gradient, 
calculating the turbulent heat and momentum fluxes over 
wide intervals of these parameters. In Section ITVl we de- 
scribe the subcritical turbulence, driven by the ion tem- 
perature gradient (ITG) and the parallel velocity gradi- 
ent (PVG), that is responsible for the turbulent fluxes. 
In Section |Vl we interpolate our results and determine 
how the temperature and flow gradients depend on the 
heat and momentum fluxes: a transport bifurcation is 
obtained as a result. In Section IVTl we construct a para- 
metric model of the transport. This allows us, in Section 
IVIIl to study the effect of the transport bifurcation on the 
temperature gradient and investigate the range of heat 
and momentum flux values for which we expect transi- 
tions to reduced transport to occur. In Section IVIIII we 
summarise and discuss the results. 



II. EQUATIONS AND NUMERICS 
A. Gyrokinetics with Velocity Shear 

We use the gyrokinetic approximation^ to calculate 
the way in which the equilibrium gradients of the tem- 
perature and flow affect the transport of heat and mo- 
mentum by turbulent fluctuations. For a detailed de- 



scription of high-flow gyrokinetics, the reader is referred 
^q 18,24,25 g^j^j references therein. Here a brief summary 
is presented. 

Wc consider a toroidal plasma with an axisymmetric 
equilibrium magnetic field B. The field lines lie on closed 
flux surfaces, which can be labelled by the poloidal mag- 
netic flux contained within each surface. The magnetic 
field can be written as J5 = Wcfi x Wtp + B^RVcj), where (j) 
is the toroidal angle, R is the radial coordinate measured 
from the central axis of the torus (the major radius), 
and Bfj, is the toroidal magnetic field. Owing to the fast 
motion of the particles along the magnetic field lines, 
equilibrium quantities such as the density and temper- 
ature are constant on each flux surface. In the current 
study, which considers the Cyclone Base Case regime^^, 
the magnetic flux surfaces are concentric and circular in 
cross-section, so that Vij] = BgRWr, where r is the mi- 
nor radius of the torus and Bg the poloidal magnetic field. 
Thus, r may be used as the flux label, and the magnetic 
field may be written a,s B ~ BgRVcj) x Vr + B^RVcj). 

We allow an equilibrium plasma flow it, of the same 
order as the ion thermal velocity 



Vthl = 




where Ti is the temperature of the ions and is their 
mass 0. It can be shown that such a flow is toroidal, as 
any poloidal component will be quickly dampe d^^'^^ . The 
flow can then be expressed as m = cuR^'S/cf), where uj is 
the angular velocity of the flow (which must be constant 
on a flux surface). 

The state of each species s of charged particles can be 
described by its distribution function, /g, whose evolu- 
tion is given by the Vlasov-Landau equation: 

where v is the velocity coordinate and C the collision 
operator. In order to describe the turbulent fluctuations 
which give rise to the transport, the Vlasov-Landau equa- 
tion can be expanded by splitting fs into an equilibrium 
part Fs and a perturbed part Sfs- The latter is further 
split into averaged and fluctuating parts; so 

fs = Fs + {Sfs) + Sfs (3) 



^ For the density, this is only true if the toroidal flow is subsonic, 
which is what we assume below, see Eq. \6\ 

^ The definition of the thermal velocity given in JTJ is chosen 
in accordance with the analytical worksilii^ cited in this paper; 
maintaining correspondence withi^iii, where the ion thermal ve- 
locity was defined as y'Ti/mi, results in various factors of in 
the normalisations of the output from simulations. 
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where the angle brackets denote a spatial and temporal 
average over fluctuations. It is assumed that the pertur- 
bations are small and that the gyrokinetic orderin g^^i^^ 
holds, so that: 



^^^^ 

F„ k \ ili L ' 



(4) 



where fc|| and k± are the typical parallel and perpendicu- 
lar wavenumbers of the fluctuations, 7 is the growth rate 
of the fluctuations, fii and pi are the gyrofrequency and 
Larmor radius of the ions, respectively, and L is scale 
length of the variation of Fg . It is then possible to show 
that to lowest order, Fs is a local Maxwellian in the 
frame of the equilibrium flow u, and that, to first order 



m e, 



Sfs--^F. 



hs {t,Rs,fis,es) 



(5) 



where ZgC is the charge of species s and (p is the perturbed 
electrostatic potential. The non-Boltzmann part of the 
distribution function, hg, is the gyrocentre distribution, 
independent of the gyrophase and a function of time t, 
the gyrocentre position Rg , the magnetic moment fig and 
the energy of the particlei^. The gyrocentre position 
Rs = r — b X v/rig, where r is the position coordinate, 
rig is the gyrofrequency of species s, the unit vector b = 
B/B is the direction of the equilibrium magnetic field, 
and B is its magnitude. 

The turbulence can be affected both by the gradient 
of the flow and by its magnitude, the latter through the 
Coriolis and centrifugal drifts. Here we neglect the ef- 
fects of the magnitude of the flow (the consequences of 
this simplification are discussed briefly in Section |VIII[) 
but allow the gradient of the flow to be of order the 
growth rate of the fluctuations, Rdu/dr ^ 7, as is nec- 
essary for the flow gradient to affect the turbulence non- 
negligibly^. 

Under these combined assumptions, it can be shown 
that hg evolves according to the following gyrokinetic 
equation: 



W|| b - 



dt 



Vd 



B^viiRduj 



dhg 
~dt ' 

ZgC 



Bn, 



dr 
1 dn, 
dr 



Tg dr \Tg 2 
+ {C[hg])j,, (7) 



where rig and Tg are the equilibrium density and tem- 
perature of species s, respectively, is a gyroaverage 
at constant Rg, d/dt = d/dt + u ■ Vjj is the mag- 
netic drift velocity, v^^ = ay^ {sg — figB) /nig is the paral- 
lel (peculiar) velocity, a is the sign of v^^ , rUg is the mass 
of species s, and c is the speed of light. We have further 
simplified the problem by assuming purely electrostatic 
perturbations {5B = 0), so the system is closed by the 
quasineutrality condition: 



E 



(8) 



where ()^ is a gyroaverage at constant r. The electrons 
are treated through a modified Boltzmann response^^: 



(9) 



where the overbar denotes averaging over the flux surface. 
Thus, Eq. ([7|) is only solved for ions: s = i. 

Knowing hi, we can calculate the turbulent heat and 
angular momentum fluxes associated with a given minor 
radius and given values of the gradients. 



Qt = 



1 c 
d^v-m^v^ — {hi 

Z D 



b X Vyj) • Vr 



(10) 



B. Numerical Setup 



(11) 



With the additional assumption that the hke-particle coUision 
frequency Vas satisfies £-^7 ^ Vaa ^ 7- 

This can be systematised by assuming that the magnitude of 
the flow is smaU but the gradient of the flow is large in a Mach- 
number ordering subsidiary to namely 



BuJ Rdw 1 

M < 1, ~ — > 1 

Vthi 1^ dr M 



(6) 



where e < M < 1. 



The nonlinear gyrokinetic equation ([7|) is solved us- 
ing the gyrokinetic code GS3^i^. Its solution depends 
on (among other quantities) the equilibrium gradients: 
drii/dr, dTi/dr and {B^/B)dLd/dr (the latter of which 
is the gradient the of parallel component of the back- 
ground flow). To simplify the treatment of these gradi- 
ents, we consider a small region of the plasma known as 
a flux tubers, which encompasses a set of magnetic field 
lines, extending to include several turbulence decorrela- 
tion lengths in each direction. This allows equilibrium 
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quantities to be expanded around a reference flux sur- 
face near the centre of the domain, labelled by rp. Thus 
we may write: 



ni (r) = Ui (ro) + (r - tq) 



riio 



drii 
dr 

1 - (r - ro) 



r=ro 
f 



Ln 



dr 



iO 



^-{r-ro) — 



duj 



cj (r) = w (ro) + (r - ro) — 

dr 



Wo + (r - ro) — 7_B. 



(12) 



(13) 



(14) 



where L„ is the equilibrium density scale length, Lt the 
equilibrium temperature scale length. The perpendicular 
velocity shear and the magnetic safety factor are: 



ro duj 

lE^ — -r- 
qo dr 



90 



RoBe 



(15) 



respectively, where Rq is the major radius of the magnetic 
axis. For a given equilibrium magnetic field, the solution 
of a flux tube simulation and, therefore, the values of Qt 
and lit, is a function of the three parameters L„, Lt and 
jE- Defining the local radial coordinate as 



x = r — ro, (16) 
and the second perpendicular coordinate y as 



BgRo 
Ba, 



(17) 



where 9 is the poloidal angl^l, then transforming to a 
frame rotating with angular frequency ujq and keeping 
terms to the correct order in the gyrokinetic expansion, 
we find for any perturbed quantity ^ : 



The gyrokinetic equation ([7]) may now be written 



^ It should be noted that because the sign of Bg = B ■ V9 in this 
paper is opposite to the sign of the poloidal field in GS2, the 
coordinate y, defined in I I17I I. also has the opposite sign to the 
GS2 coordinate of the same name. 



dhi dhi ( 



dt 



dt 



9{v)r^ ^ B^v\\Rq^^E d{^)R^ 
dy 



B dy 



Lr. 



Bil, 



ro dy 
3^ 



F,, 



In this investigation, we vary the perpendicular flow 
shear je, and k = Rq/Lt- Other parameters are kept 
fixed: 



Ln 



2.2, qo = 1.4, 



i?o 



0.18. 



(20) 



The density gradient, the magnetic safety factor and the 
aspect ratio are chosen to conform to the Cyclone Base 
Case, foUowingiS andii, while, as we explained above, 
the magnetic shear is set to 0. 

Note that the value of qo effectively controls the 
strength of the PVG drive (the term proportional to qojE 
on the right hand side of for a given velocity shear 
7b. The relatively low value of qo in the Cyclone Base 
Case reduces the destabilising effect of the PVG — com- 
pared for example to the Waltz Standard Case where 
go = 2.02. 

Collisions are included by means of a model collision 
operator, which includes the scattering of particles in 
both pitch angle and energy^Si^. The numerical ion-ion 
collision frequency i/^ = O.Olvthi/ Ro- 

The effects of the perpendicular flow shear are in- 
cluded in the code as follows^. Expanding ^ = 
J2k k ^i'^jt) exp[i {kyy + kxx)], where z is the coordi- 
nate along the field line, the perpendicular flow shear 
terms in (jl9p can be eliminated by adding time variation 
to the radial wavenumber: 



kxit) = kxo - jEkyt, 



so that 



^+X^E — 

dt dy 



^ i[k^{t)x+kyy] 

dt 



(21) 



(22) 



All fluxes will be reported in dimensionless units by 
normalising them to the gyro-Bohm values 



QgB 



2^/2i?g ' 



n„ 



4i?o 



(23) 



The flow shear will be everywhere normalised to vthi/Ro] 
hence at this point we rescale je accordingly: 



7-E 



Ro 



IE- 



(24) 
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The resolution of the majority of simulations was 
64 X 32 X 14 X 24 X 8 X 2 — the number of gridpoints in the 
kxo and ky spectral directions, in the z spatial direction, 
the average number of pitch angles (i.e. the number of 
values of for a given , which varies with the parallel 
coordinate), the number of gridpoints in energy space Si 
and the two values of a. This grid provided sufficient 
scale separation in both spatial directions perpendicular 
to the magnetic field for calculating the turbulent trans- 
port, and sufficient resolution in velocity space to calcu- 
late the velocity integrals in Maxwell's equations with the 
required accuracy. A discussion of the parallel resolution 
is given in Section IIVI 



III. TURBULENT FLUXES 




A. Heat Flux 



FIG. 1. ITG instability and turbulence at zero flow shear: (a) 
linear growth rate; (b) turbulent heat flux vs the normalised 
temperature gradient k = Ro/Lt- The case of k = 11, further 
explored in Fig. |3] is marked by *. 



Using simulations in the manner described in Section 
im the heat and momentum fluxes were calculated for k 
values in the range 4 < k < 13, vs flow shear values in 
the range < 7^ < 2. 

Fig. [IJb) shows the heat flux vs the ITG in the absence 
of flow shear. When = 0, the critical temperature 
gradient above which there is ITG-driven turbulence is 
K — 4.4. When k is increased above this threshold, the 
heat flux increases rapidly up to more than a hundred 
times the gyro-Bohm value. 

Fig. [5] shows the turbulent heat and momentum fluxes 
vs the flow shear for different ITG values. As the flow 
shear is increased from 0, the heat flux initially responds 
weakly, either increasing or decreasing slightly. As the 
flow shear je approaches 1, the heat flux dips sharply 
for all values of the ITG. For moderate temperature gra- 
dients {k ^ 11), the turbulence is suppressed altogether 
and the heat flux drops to zero. The suppression of tur- 
bulence also happens at finite magnetic shearsiS^ but for 
a significantly narrower range of k H. 

For larger je, the heat flux starts to rise again. This 
increase is due to the PVG — we have verified that this 
effect disappears if the term proportional to qq^e on the 
right hand side of eq (fT9|l is artificially set to 0. This 
revival of the turbulent transport at large shears was 
also observed ini^, but was absent from the quasilinear 
study conducted in^, which showed the turbulence being 
completely suppressed above a sufficiently large toroidal 
shear. 
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FIG. 2. Turbulent heat (a) and toroidal angular momentum 
(b) fluxes (normalised to gyro-Bohm values) as functions of 
flow shear for different values of the ion temperature gradient 
k; (c) turbulent Prandtl number vs flow shear (for cases where 
the heat flux is non-zero). 



B. Momentum Flux 



This is partly owing to the fact that at finite magnetic shear 
growing linear eigenmodes exist for non-zero flow shear, whereas 
at zero magnetic shear there are no such eigenmodes except when 
the flow shear is also zero; this is discussed further in Section llVI 



The momentum flux is zero when 7b = (as might 
be expected in an up-down symmetric case,— i^). As 
7e increases, so does the momentum flux — a result of 
the turbulent viscosity. However, when the flow shear 
reaches values at which it starts to suppress turbulence 



6 



significantly, the trend is reversed and a remarkable sit- 
uation arises in which increasing flow shear reduces the 
transport of momentum. This behaviour persists until 7_e 
reaches even larger values, and the turbulence is reignited 
by the PVG drive. The positive correlation between the 
velocity shear and the momentum flux is then reestab- 
lished. It is the existence of the window of suppressed 
momentum transport at moderate and k that will 
enable the transport bifurcation analysed in Section IVl 



C. Turbulent Prandtl Number 

There is a clear correlation between the heat and the 
momentum flux, which is best quantified in terms of the 
turbulent Prandtl number 



Prt = — 

Xt 



Qt/QgB IE V^RoQo 



(25) 



where the turbulent viscosity and the turbulent heat dif- 
fusivity are 



Xt 



lit 1 VthiroPi 
Qt 1 vthipf 

QgB K 2V2R0 ' 



(26) 
(27) 



respectively. 

The Prandtl number obtained from our simulations is 
plotted in Fig. [D^c). There is a similar basic trend for all 
values of k: Pit rises from approximately 0.5 when je — 
0.1 Q, peaks at 7^; ~ 1 and then drops to approximately 
0.6 when 7^; ~ 2. For those intermediate values of 7^ 
where the turbulent fluxes arc reduced or tend to 0, Prt 
rises sharply, reaching just above 0.7 for low values of k. 
The location of this sharp rise varies with temperature 
gradient similarly to the location where the turbulence is 
suppressed. 

Although the Prandtl number does vary with both je 
and K, this dependence is relatively weak compared to 
the individual dependence of vt and Xt on these param- 
eters. Thus, approximating Prt ^ 0.55 everywhere keeps 
it within approximately 25% of the true value for the en- 
tire range of 7^. This will prove convenient in construct- 
ing a model of turbulent transport presented in Section 

m 



^ The low value of Prt for small 7^ has been observed before 
and is sometimes referred to as a shear pinch^^i^. It occurs 
because, at low 7^, the perpendicular flow shear can give rise to 
a contribution to the viscous stress that has the opposite sign to 
that of ITG-driven turbulence with zero flow shear, reducing the 
overall diffusive transport. 



IV. SUBCRITICAL TURBULENCE 

Before studying the implications of the results pre- 
sented above, we first discuss the nature of the fluctu- 
ations that give rise to the behaviour of the turbulent 
fluxes discussed in Section Hill 
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FIG. 3. Transiently growing linear modes a.t k = 11. (a) The 
heat flux as a function of time, normalised to its value at to, 
where to — 0.06\/^Ro/vthi; so chosen to skip the short initial 
transient associated with a particular choice of initial condi- 
tion. All modes grow transiently for 7_b > 0. (b) Duration 
of transient growth r-y from t = to to the peak value of Qt vs 
flow shear. 

With jE — 0, the ITG instability exhibits robust linear 
growth, with a threshold of k = 4.4 (Fig. [H^a)). However, 
for any finite value of the flow shear there are no growing 
linear eigenmodes in the system: all modes grow only 
transiently before decaying (Fig [31(a)). While formally 
this means that 7^; — > is a singular limit, there is is in 
fact no physical discontinuity in behavior: as 7_e — >■ 0, 
the duration of the transient growth, Tj, tends to infinity 
(roughly as 1/7_e; see Fig. |31^b)). There is therefore a 
continuous transition at 7^; = from a transient mode 
that grows for an infinitely large time, to a growing eigen- 
mode. The lack of growing eigenmodes when 7^; > 
stems from the wavenumber drift (|2ip : as time tends 
to infinity at finite flow shear and zero magnetic shear, 
the radial wavenumbers increase inexorably through the 
shearing of the modes, until they are extinguished by 
finite-Larmor-radius effects. 

This situation differs somewhat from the case with fi- 
nite magnetic shear S^^^, where there are linearly grow- 
ing modes for a finite range of je- At finite magnetic 
shear, growing eigenmodes are possible because while je 
leads to an effective dependence of the radial wavenum- 
bers on time, s makes them dependent also on the po- 
sition along the field line. Therefore, for modes moving 
with a speed proportional to Je/s, the magnetic shear 
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FIG. 4. Two measures of the strength of the transiently grow- 
ing linear modes at k = 11: (a) the effective growth rate of 
the mode and (b) the number of exponentiations of the heat 
flux during the growing phase, both vs. the flow shear. The 
apparent discontinuity at 7_b = in (a) is a result of taking 
the average growth rate of the heat flux over < t < r.y 
rather than the initial or peak growth rate; for all 7_e > the 
average will include a period where the growth rate tends to 
zero. 



cancels the effect of the velocity shear on kx (t)^^. When 
the magnetic shear becomes very small, or the flow shear 
very large, the required mode velocity becomes much 
greater than the thermal speed of the ions and the can- 
cellation is no longer possible because the mode cannot 
travel so fastis. As a consequence, growing linear eigen- 
modes only exist for non-zero magnetic shear and only 
up to a finite value of the flow shear—. 

In a standard picture of ITG turbulence with 'je = 
the turbulence is driven by unstable linear modes. With 
the exception of a narrow interval of temperature gradi- 
ents where self-generated zonal flows suppress the turbu- 
lence (the Dimits shift), the presence and the amplitude 
of the turbulence are largely determined by the presence 
and the growth rate of those unstable modes. In the 
present case, by contrast, we have strong levels of turbu- 
lence sustained nonlinearly in a parameter regime where 
there are no linearly unstable modes, a phenomenon first 
discovered in tokamak turbulence im^ andii, but well 
known in various hydrodynamical contexts^. This phe- 
nomenon is known as subcritical turbulence. 

Subcritical turbulence differs from standard 
instability-driven turbulence in several important 
ways. Firstly, because there is no linear instability, 
turbulence will not grow from initial perturbations of 
arbitrarily small amplitudoi^. Fluctuations must be 
initialised with sufhcient amplitude (generally of the 
order of their amplitude in the saturated state) in order 
for turbulence to be sustained; thus, the absence of 
sustained turbulence in a numerical experiment may 
merely indicate that the initial amplitude is insufficient. 



not that the plasma is quiescent. 

The second difference concerns the scales at which the 
turbulent energy resides. In ITG-driven supercritical tur- 
bulence (7^; = 0), these scales are those where the lin- 
ear drive injects energy — this tends to correspond to 
/Cii qoRq ^ 1 and kyPi relatively narrowly concentrated 
around values of a fraction of unity (at least for low val- 
ues of qo and moderate temperature gradients; see^ and 
Fig. El^a)). In subcritical turbulence, the preferred wave- 
lengths appear to be those at which the amplification of 
the transient modes is strongest. In the case of turbu- 
lence where the PVG drive is dominant, Ref has shown 
that the amplification is maximised along a line in Fourier 
space where 



kyPi 



ki\Ro 



(28) 



Fig. [SJb) shows the spectrum of strongly PVG-driven 
turbulence at high fiow shear (7^ = 2.2). We sec that 
although the result (pS)) is based on linear theory and is 
derived for slab geometry, it appears to describe the peak 
of the spectrum reasonably well. The spectrum extends 
to much higher parallel wavenumbers than in the stan- 
dard ITG case; consequently, higher parallel resolution is 
necessary to resolve i1o. 

Finally, faced with subcritical turbulence, we are left 
without an intuitively obvious way of estimating its satu- 
ration level and, consequently, the turbulent fluxes. It is 
the presence of linear eigenmodes with a defined growth 
rate (positive or negative) which has enabled the quasilin- 
ear modelling of the heat flux in many situations without 
resorting to full nonlinear simulations (— was an early ex- 
position; see^ for an overview and^i^ for recent work) . 
When the growth of all modes is transient, such models 
will not work in their current form. The question arises 
as to which characteristics of the linear transient growth 
are relevant in the resulting nonlinear state. The prob- 
lem is as yet unsolved, but here we consider two natural 
measures of subcritical driving. 

Let us flrst dcflne the effective transient growth rate^ 



7eff 



i.^^g.(i 



2t^ Qt{t = to)' 



(29) 



This effective growth rate is plotted in Fig. IHa). It 
initially decreases with increasing flow shear, but then 



* Studies of the effect of increasing the parallel resolution showed 
that at values of flow shear 7^ ^ 1-2, the parallel resolution used 
for the parameter scan described in Section IlIII (14 gridpoints) 
was insufficient and led to errors in the critical temperature gra- 
dient (above which subcritical turbulence could be sustained) of 
order 5-10%. We do not consider this error significant enough to 
merit a repeat of the parameter scan with substantially higher 
parallel resolution (although in principle such an exercise would 
be useful) . 
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starts to increase for 7^ > 0.4 as the PVG drive becomes 
significant. The effective growth rate seems to correlate 
with the behaviour of the heat flux (Fig. El^a)): it first 
decreases and then increases with increasing flow shear. 
However, the minimum in the effective growth rate (for 
K = 11 as shown in Fig. HJ^a)) occurs at 7^; ~ 0.4 whereas 
the minimum in the heat flux for the same k occurs at 
7-E ~ 1- 

Besides the rate of growth, what is likely to be deci- 
sive in determining whether the finite transient growth of 
perturbations proves sufficient to maintain a steady level 
of turbulence is the time that this growth lasts. A nat- 
ural diagnostic that accounts for both effects is the total 
transient amplification of the modes before they start to 
decayii^di: 

N = r^7eff. (30) 

It drops precipitously as 7^ is increased from 0, reaches a 
minimum at 7^ = 0.6 (which does not coincide with ~ 
1 , where Qt is minimal) and then appears to increase very 
gently with 7_e El — in sharp contrast with the very rapid 
increase of the nonlinear fiuxes at high fiow. 

Thus, while the linear behaviour shows that the non- 
linear fiuxes are somewhat determined by the vigour of 
the underlying linear drive, a quantitative model of sub- 
critical PVG turbulence that would allow the heat fiuxes 
to be predicted from the characteristics of transient linear 
growth remains an unsolved problem. 



1.6 rn 1 




-10 -5 5 

fc||(30-Ro 



FIG. 5. The spectrum of turbulent fluctuations for nor- 
mal and subcritical turbulence: fcy X]fe I'^fcxo.fcyl^ (arbitrary 
units) vs ky and for (a) 7_e ~ 0.0, k = 12 (ITG turbulence 
with no flow shear) and (b) je = 2.2, k = 12 (subcritical, 
strong PVG-driven turbulence)0. Also shown in (b) is the 
line of maximum transient amplification H28p . as calculated 
iftii. 

^ It should be noted that as the sign of ky is opposite in GS2 to 
the present work, these graphs were plotted via a parity 
transformation. 



^ This behaviour approximately agrees with the analytical result 
o£i4, which predicts N — >■ 0.225goflo/ro ~ 1.75 as 'je — >■ 00. 



V. TRANSPORT BIFURCATION 
A. Possibility of Bifurcation 

In Section ITTTl we demonstrated the existence of a wide 
interval in 7^; in which a sheared toroidal equilibrium fiow 
completely suppresses turbulent transport. If such a sup- 
pression could be achieved experimentally in a tokamak, 
confinement of energy would be dramatically improved. 
Unfortunately, while it is possible to specify the fiow 
shear in numerical simulations, in experiment it can only 
be varied indirectly by applying torque, and the effect of 
that torque is strongly dependent on how quickly the an- 
gular momentum escapes from the plasma. If only a lim- 
ited amount of torque can be injected and the momentum 
flux is too large, it could be impossible to achieve fiow 
shears that are large enough to suppress the turbulence. 
FigHfb), however, suggests an intriguing possibility. For 
all temperature gradients, the momentum flux at flrst 
increases, reaches a maximum and then decreases before 
increasing again at higher fiow shears. There is, there- 
fore, a parameter window in which increasing the fiow 
shear decreases the transport of momentum. If this were 
to happen, the momentum would build up, increasing 
the fiow shear, which would further decrease the trans- 
port of momentum. Such a positive feedback could lead 
to a bifurcation and so it might seem that high values 
of fiow shear could be reached without excessive input 
of momentum, a possibility which was discussed in^Fl. 
However, Fig [2Kb) shows the momentum fiux at constant 
ion temperature gradient, which would not, in fact, stay 
constant during this process. Indeed, as soon as fiow 
shear began to suppress the turbulent transport of mo- 
mentum, it would also suppress the turbulent transport 
of heat, causing an increase in the temperature gradi- 
ent as well as in the flow gradient. This increase in the 
temperature gradient would restore the turbulence to its 
former levels. Nevertheless, we will show in this section 
that when neoclassical transport is also taken into ac- 
count, a bifurcation is possible. 



B. Inverting the Problem 

What can actually be controlled in a steady-state ex- 
perimental situation is the flow of heat and momentum 



Note that although a similar mechanism was discussed ina», that 
study considered flow shears up to a maximum of approximately 
0.08cs/a, where Cs = ^2Telmi is the sound speed and a is the 
minor radius of the plasma (in this study Cs = vthi and a = -Ro). 
Thus they were considering a maximum in the momentum flux 
which occurs at low flow shear and finite magnetic shear, and 
which is discussed in more detail irki£. As can be seen from 
Refpifi, the corresponding jump in flow shear is much smaller, and 
the range of k values at which such a maximum in the momentum 
flux exists is much narrower, than in the present case (Fig. [2{b)). 
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through a particular surface. This means that we must 
switch from using 7^ and k as independent parameters 
to using the total fluxes of heat and momentum, Q and 
n respectively. 

Let us consider an experimental set-up where the heat 
and the momentum are injected by beams of neutral par- 
ticles. We assume that the energy and momentum from 
those neutral beams are deposited uniformly across the 
torus. We further assume that the beams are tangential 
to the magnetic axis and we ignore corrections of order 
the aspect ratio of the device. Then: 



1 VfNbm,vl/2 2^/2i?oPNBI 



QgB QgB 47r2roi?o 

n 1 VfNi,m,Vb 



An^roniTiVthiPi 



n/n^ 



n 



gB 



Vtht 



Q/QgB V^Vb V^^'NBI/ 



1/2 



(32) 



(33) 



where Nb is the number of neutral beam particles in- 
jected per unit time, Vb is the beam particle velocity, 
^^NBi = Nbmivl/2 is the beam power, Snbi — 'niivl/2 is 
the beam particle energy and V/ is the volume fraction of 
the plasma enclosed by the flux surface. Thus, the total 
heat flux Q is determined by the beam power Pnbi ([ST]) , 
whereas the ratio of the total momentum angular flux 
to the total heat flux, II/Q, is determined by the beam 
particle energy i?NBi ((33)) . We want to know whether, by 
varying Q and II/Q, it is possible to reach, or to trigger 
a transition to, a high-flow-shear regime where the tur- 
bulent transport is suppressed. Thus, it is necessary to 
convert the dependence of Q and 11 on k and 7^, deter- 
mined from local simulations, to a dependence of n and 
7_B on Q and II/Q. 



C. Interpolation 

The gyrokinetic code GS2 gives the fluxes as a function 
of K and 7b. Given the expense of the nonlinear simula- 
tions necessary to do this, it is computationally challeng- 
ing to use GS2 as a root finder to invert the problem and 
find the gradients as a function of the fluxes^i^. Instead 
we interpolate within the set of data points described in 
Section IIIII (as well as additional simulations in param- 
eter regions of particular interest) to obtain the fluxes 
as continuous functions of the gradients; these functions 
can then be inverted numerically to give k and as 
functions of the fluxes. 

Interpolation in a multidimensional parameter space 
is a nontrivial operation. A standard technique is to use 
radial basis functions^, which weigh each data point by 
its distance in parameter space from the point of interest 
(after the parameter space has been normalised to ensure 
that variation occurs on the same scale in each coordi- 
nate). There are many choices of the function, or kernel, 
which is used to calculate the relative importance of each 



data point. Here we choose a linear kernel (equivalent to 
linear interpolation in the case of only two points)^. Us- 
ing this, the values of the fluxes at a point p — {-fE, k) 
can be calculated as follows: 



(p) = E 



(34) 
(35) 



where pj arc the input gradients for the nonlinear simula- 
tion labelled j, and the weights w'^ and arc calculated 
so as to satisfy for all j: 

Ot(Pj) = Otj ,nt(pj) = n,j, (36) 



where Qtj and Iltj arc the values of the fluxes obtained 
from simulation j . 



D. Neoclassical Transport 

If the turbulent transport is successfully suppressed, 
neoclassical (coUisional) transport becomes important. 
The total fluxes arc the sum of the neoclassical and tur- 
bulent contributions: 



Q = Qt + Qn, n = nt + n„, 

where the neoclassical fluxes are calculated as 

Qn _ 2V2R0 
QgB VthiPi 

n„ ARlqo 

- VnlE- 



(37) 

(38) 
(39) 



IlgB roVtMPi 

and the neoclassical thermal diffusivity and viscosity are 



R 



Xn ^ 0.66 ( -- ) q^pfiy.. 



3/2 



(40) 

~ O.lqlp^^Vu. (41) 



The formulae for the neoclassical diffusivity and viscosity 
in the large-aspect-ratio, concentric-circular-flux-surface, 
banana (lya ^ vthi/Rq) regime, were taken froroii. Iitii, 
the ion-ion collision frequency is defined as: 



/2^Zfn, e^ In A 



(42) 



(where In A is the Coulomb logarithm). Thus, the neo- 
classical transport is a function of the temperature and 
density, and scales with these quantities in a different 
way to the turbulent transport. To determine the ra- 
tio between the neoclassical and turbulent transport it 
is necessary to choose specific values for the temperature 
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and density. Here, as inii, we take i^u = ly^ /2. Later, 
when we consider the case of specific tokamaks, we will 
use typical values of Ui and Ti to estimate i^u. 

Using the results PO]) and (PT|) the neoclassical Prandtl 
number can be calculated from the parameters given in 
(PU]) and is found to be 



Pr„ = — ~ 0.01 < Pit. 

Xn 



(43) 



Thus the neoclassical Prandtl number is much smaller 
than the turbulent Prandtl number (see Fig. [DJc)). In 
other words, the neoclassical transport of momentum is 
much smaller than the neoclassical transport of heat, in 
contrast to turbulent transport, for which the fluxes of 
momentum and heat are comparable. While the formulae 
(HI]) and (HOI) are approximate, we emphasise that the 
qualitative result of the following section is not affected 
by small changes in the values of Xn and ;/„, provided that 
the property Pr„ ^ Pr^ continues to hold. In particular, 
this means that this qualitative result is not affected by 
small changes in the value of va . 



E. Bifurcation 
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FIG. 6. Momentum flux divided by the heat flux (each nor- 
malised by the corresponding gyro-Bohm estimate) vs the flow 
gradient at a constant value of the heat flux Q = 2.6QgB 
plotted using (a) interpolation from the data explained in 
Section IV CI and (b) the parameterisation from Section IVII 
Also shown are the neoclassical contribution to the momen- 
tum flux (dotted line) and the momentum flux at constant 
heat flux without the neoclassical contribution (dashed line). 



In Fig. ISl^a), the momentum flux is once more plotted 
against ^e, but this time keeping the heat flux constant 
at Q/QgB ~ 2.6. The local maximum in the momentum 
flux still exists. Starting at this local maximum (point 
A), if the torque on the plasma was to be increased at 
constant Q, the flux of momentum could only increase 
if the flow shear were to jump to a much higher value 



(point B) where the PVG instability would drive turbu- 
lent momentum flux. A bifurcation is manifest. 

The inclusion of the neoclassical fluxes is critical in 
obtaining this result. To illustrate this. Fig. |6] also 
shows the momentum flux at constant heat flux with- 
out the neoclassical contribution to the fluxes, i.e., with 
Xn = Vn = 0. If, in such a synthetic "purely turbulent" 
transport case, the flow shear is increased at constant Q, 
the temperature gradient increases to maintain the tur- 
bulent heat flux, and, as shown in Fig. [6l II/Q increases 
monotonically with je- The key property of neoclassical 
transport that helps change this behaviour and produce 
a bifurcation is that the neoclassical Prandtl number is 
much smaller than the turbulent Prandtl number. This 
means that as turbulent transport is suppressed, the sys- 
tem goes through a regime where the neoclassical contri- 
bution to the heat flux is significant, while the neoclas- 
sical contribution to the momentum fiux is not. So as 
n/Q is increased at constant Q, the turbulence is sup- 
pressed, and the temperature gradient starts to rise, but 
as this happens more of the heat fiux is transported neo- 
classically, and so the feedback loop breaks down: it is 
no longer necessary to increase the levels of turbulence 
to maintain the same Q. The same is not true of the mo- 
mentum: the neoclassical viscosity cannot make up for 
the lack of turbulence, and so the transport of momen- 
tum peaks and then falls. 

At high flow shears the turbulent momentum flux in- 
creases rapidly once again; this turbulent flux is driven 
by the PVG, as discussed in Section Hill As a result, we 
observe that the bifurcation results in a turbulent state 
at B, with a signiflcant turbulent momentum flux. This 
is a substantial difference from the situation envisioned 
ii*^, where a transport bifurcation was predicted using a 
reduced quasilinear model. Their model did not (and, 
being a quasilinear model, could not) predict the exis- 
tence of the PVG-driven subcritical turbulence at high 
flow shears; instead it predicted a bifurcation resulting 
in a non-turbulent state where all transport was neo- 
classical. Thus a full nonlinear analysis is necessary to 
describe accurately the reduced transport state produced 
by the bifurcation that we find in a turbulent plasma. 

Restoring dimensions, we find that the bifurcation 
results in a jump in the flow shear from 0.25 vthi/Ro 
to \.n Vthi/ Rq- It is instructive to consider whether 
such values of shear appear in current devices. Taking 
the measured values of the perpendicular flow shear, Ti 
and i?o, in high-performance discharges in the JET and 
MAST tokamaks^Sii^, we estimate that flow shears of up 
to approximately Vthi / Ro have been recorded in both de- 
vices, and thus that values of shear of the order examined 
in this paper have been observed. 



VI. PARAMETERISED MODEL 

Through simple interpolation, without making any as- 
sumptions about the way the fluxes depend on the gradi- 
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ents, we have shown the existence of a bifurcation in the 
gradients at constant fluxes. However, even interpolat- 
ing one line of constant Q requires a very large number 
of data points in the region of the line. To produce the 
required data set for Fig. [6l^a), we had to perform about 
350 nonlinear gyrokinetic simulations, each run to satu- 
ration. In order to extend our understanding of the bifur- 
cation, and of the range of flux values for which similar bi- 
furcations can occur, we will first consider the behaviour 
of the turbulent fluxes in further detail and then use this 
analysis to construct a simple parameterised model of the 
turbulent fluxes Qt and lit as functions of and k. 



A. Modelling Qt 

1. Dependence of Qt on k 
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FIG. 7. Turbulent heat flux vs the temperature gradient for 
diiTerent values of the flow shear, showing (a) low-shear values 
7b < 1, (b) a close up of the low Qt region in (a), and (c) 
high-shear values 7b > 1 

We wish to construct a parameterised model of Qt as 
a function of and k. In Section IlIII we described the 
dependence of Qt on 7^; we now consider the dependence 
of Qt on K. Both experimentallji^ and numericalljsi^, it 
is usually found that Qt increases very sharply with k — a 
property known as stiff transport. A recent experimental 
study^ has indicated that increasing the flow shear at 
low magnetic shear might have the effect of reducing the 
stiffness. Fig. [7] shows that the effects of flow shear on 



dQt/dn are in fact quite complex. Let us consider the 
cases of low flow shear, 7^ < 1, and high flow shear, 
7b ?i Ij separately. 

For the case of 7^; < 1, shown in Fig. El^a), the thresh- 
old in K above which turbulence can be sustained non- 
linearly increases rapidly with 7^, as the perpendicular 
shear suppresses the ITG instability. Above this thresh- 
old there are broadly three regions: low, intermediate 
and high Qt- 

At high Qt, far above the threshold, for any value of 
flow shear the heat flux eventually asymptotes to the 
same dependence as it has at ^e — 0. In other words, 
as the ITG drive becomes very strong, the effect of flow 
shear becomes negligible. 

At intermediate Qt, the heat flux rises rapidly from 
low values to join the universal, high-Qt asymptotic. As 
7b increases, because the threshold rises, the heat flux 
rises more rapidly with k above the threshold; thus at 
intermediate values of Qt, the stiffness dQt/dn increases 
with jE', this is shown in Fig. [S] 

At low Qt, just above the threshold, the heat flux rises 
very slowly, that is, the stiffness dQt/dn is very low (see 
Fig. |8]). This low Qt region is shown in detail in Fig. 
[Tt^b). Such is the sharpness of the transition from the low- 
stiffness low-Qt region to the high-stiffness intermediate- 
Qt region that there appear to be two distinct thresholds. 
The first threshold is the transition from no turbulence to 
non-zero turbulent transport. Above the first threshold 
turbulence is present but Qt rises slowly with k (the low- 
Qt region). Above the second threshold Qt rises rapidly 
(the intermediate-Qf region). These thresholds are plot- 
ted in Fig. ini The low-stifi'ness region only exists for 
< 7e < 0.8: between 0.4 < 7^ < 0.8 the first threshold 
joins the second and the low-stiffness region disappears. 

At high flow shear, 7^ > 1, there are two principal 
differences to the case with low fiow shear. Firstly, the 
low-stiffness region is not present, and there is only one 
threshold. Secondly, the critical temperature gradient for 
turbulent heat transport starts to decrease as the PVG 
drive rc-enforces the ITG drive. When je = 1-8, the 
PVG drive is strong enough to drive turbulence at very 
low k: the threshold drops to zero. As the threshold 
drops, the transport stiffness at intermediate values of 
Qt decreases, in a mirror image of the case at low flow 
shear. At high Qt, the heat flux still asymptotes to the 
universal je ~ curve. 

Thus, increasing je can both increase and decrease the 
nonlinear thresholds, and both increase and decrease the 
stiffness. The rise and fall both in the thresholds and 
the stiffness can also be seen in the finite-magnetic-shear 
simulations ofi^. The principal differences between the 
zero-magnetic-shear case considered here and that case 
are that we have a higher value of the critical gradients 
for all jE, and that we find a low-stiffness region at low 
values of Qt — a feature that seems to be absent at finite 
magnetic shear. 



12 



(5 



60 

50 
40 
30 
20 
10 




Intermediate Qt 



Low Qt 
* * f 



0.2 0.4 0.6 0.8 1 

7E 



1.2 1.4 1.6 



FIG. 8. The dependence of the heat transport stiffness dQ/dn 
on the flow shear in the \ow-Qt and intermediate-Qt regions, 
measured using simulations close to both thresholds (points, 
see Fig. [7]) and using the parameterised model of Section 
IVI A 21 (lines'). 
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FIG. 9. The first and second critical thresholds, Kci and Kc2, 
measured using simulations close to both thresholds (points, 
see Fig. [7]), and using the parameterised model of Section 
IVI A 21 (lines'). 



2. Parameterisation of Qt 

In^^ a simple model was used to characterise the be- 
haviour of the turbulent heat flux, and describes the qual- 
itative behaviour of the heat flux reasonably well. Here, 
we describe a more complex model that reproduces most 
of the features described above and is quantitatively close 
to the interpolated fluxes. We consider only the low and 
intermediate Qt regions, as the bifurcation occurs near 
the boundary between these regions. In order to de- 
scribe Qt in these regions, we have to parameterise the 
two thresholds and the transport stiffness dQt/dn. 

We parameterise the first and second critical thresholds 
as linear and quadratic functions of je respectively: 



K-cl 



Kf) + ai7_E, 



(44) 
(45) 



where kq = 4.4 is the nonlinear threshold for turbulence 
at 7b = and the parameters ai — 10.1, a2 — 17 A and 
Q!3 = —17.0 are chosen to fit the data. The modelled 
thresholds are plotted next to the measured thresholds 
in Fig. |9l As with the observed thresholds, the first 
threshold joins the second between 0.4<7£;<0.8. 

Next we parameterise the transport stiffness dQt/dti. 
The measured values of dQt/dn in the first and second re- 
gions are shown in Fig. [51 Between the first and second 



thresholds (the low Qt region). d{Qt/ Qqb) / is mod- 
elled as a constant, a^. In the intermediate Qt region, 
remembering the observation that dQt/dK broadly rises 
and falls with the nonlinear thresholds, we allow the stiff- 
ness to depend on Kc2- Thus, 



1 dQt 

QgB dK 

1 dQt 

QgB dK 



a4 Kcl < K < Kc2, (46) 

a5 + aeKc2 k > Kc2- (47) 

8.0. 



The parameters are = 1.5, ~ 3.0 and ag 
The model of dQt/dn is shown in Fig. [51 

Thus Qt is parameterised as a piecewise linear function 
of K. It is zero below the first threshold, has gradient 04 
above the first threshold and gradient as 4- aeKc2 above 
the second: 



Qt 



a4 {k - Kcl) 



K < Kcl, 

Kcl < K < Kc2, 



QgB 



Cti {Kc2 — Kcl) 

+ {a^ + aeKc2) {k - Kc2) K>Kc2- 



(48) 



The model is compared with the data points in Fig. 1101 
To reproduce the bifurcation of Section IV El in a quanti- 
tatively correct manner, it is in fact sufficient to represent 
accurately the two thresholds and the low stiffness region. 
However, by also parameterising the variation of dQ/dK, 
we have provided a model which uses six parameters to 
describe the (complicated) behaviour of the heat flux over 
a wider parameter regime with reasonable accuracy. 
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FIG. 10. The modelled heat fiux (lines) shown along with the 
simulated data points for (a) low flow shear and (b) high flow 
shear. Legends as in Figures [7j a) andlTJc). 
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B. Modelling Ut 

While the turbulent Prandtl number does vary with 
Stiid K, as discussed in Section IIII C[ this variation 
is relatively weak. Thus we model Ut by assuming a 
constant turbulent Prandtl number: 



Pit = 0.55. 



(49) 



IlgS QgB K ro 

The same choice was made in^^. 

C. The Modelled Bifurcation 

In Fig. injb) the model is used to replot the angular mo- 
mentum flux at constant Q/Qqb ~ 2.6. It has the same 
shape as the interpolated curve in Fig. IHJa), and is quan- 
titatively very close to it at low flow shear < 1.0). 
The agreement at higher flow shears is less good, which 
reflects the fact that the second critical gradient is not re- 
ally a quadratic. Nonetheless, we consider this degree of 
precision adequate. The model is used below to explore 
further properties of the transition without the need for 
extremely large numbers of additional simulations. 



VII. THE REDUCED TRANSPORT STATE 

In Section|Vl a transition was described leading to a re- 
duced transport state. Here we use the parametric model 
that was designed in Section lVTl to describe the properties 
of this state. 
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A. Heat Flux at Constant I\/Q 




Fig. [S] showed the effect of varying 11/ Q whilst keeping 
Q constant. In contrast, Fig. [TT] demonstrates the effect 
of varying Q whilst keeping 11/ Q constant. At high Q 
and low temperature gradient (marked (I) in Fig. Illf b)). 
increasing the heat flux has the effect of increasing the 
temperature gradient as might be expected: this is the 
ordinary turbulent regime. However, as Q/Qgs is low- 
ered below ^ 5 there arises a counterintuitive situation 
where decreasing Q has the effect of raising the temper- 
ature gradient. This is of course the combined effect of 
the neoclassical Prandtl number and the flow shear as de- 
scribed in Section |Vl and also discussed in detail in2^. In 
this region (marked (II) in Fig. [TlT b)). Qt becomes com- 
parable to Qm but as the heat flux is lowered, the system 
cannot drop to a neoclassical state because of the small- 
ness of the neoclassical transport of momentum. Thus 
the constancy of the ratio H/Q (i.e., the large finite flux 
of momentum) keeps the system in a turbulent state; the 
temperature gradient and the flow shear rise, and the 
curve representing Q vs k becomes flatter and curls up 



FIG. 11. Heat flux Q vs the temperature gradient k at a 
constant ratio of the momentum flux to the heat flux n/Q 
(in units of \^RQ/vthi) plotted using (a) interpolation from 
the data fSec IV Cll and (b) the parameterised model (Sec lVIl) . 
Also shown is the neoclassical contribution to the heat flux. 
Interpolation in (a) is impossible near this neoclassical line, 
where the contours are closely spaced and H/Q is multivalued, 
(b) also shows the maximum possible temperature gradient at 
low Q for a given Tl/Q (labelled "max k"). Fig. (c) plots the 
same curves as (b), showing both Q against 7_b and k, and the 
curves projected on the "/e-k. plane, illustrating the increase 
in the flow shear along each curve of constant H/Q. Points A 
and B in (a) correspond to points A and B in Fig. Ha). 



in order to maintain its distance from the neoclassical 
line. At large temperature gradients and flow shears (re- 
gion (HI) in Fig. ITlT b)). the momentum flux increases 
rapidly due to the PVG drive and the high flow gradient 
(see Fig [D^b)) and so the curve rolls over and resumes 
its original downward trend. Eventually the heat flux 
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asymptotes to the neoclassical value. 

The last part of this trend, while physically obvious, 
can only be seen using the modelled heat flux (Fig[TTJb)), 
as it is not feasible to interpolate in the region near the 
neoclassical line where the contours of constant li/Q are 
very closely spaced and II/Q becomes multivalued. 



that the turbulence levels in the reduced transport state 
are small. The dominance of 11 j over n„ given the low 
levels of turbulence reflects the fact that the flow gra- 
dient is large and the neoclassical momentum transport 
is very low compared to the neoclassical heat transport 
(Pr„ < 1). 



B. Temperature Gradient Jump 



D. Bifurcation by lowering Q/QgB 



Interpolation cannot yield the temperature gradient af- 
ter the transition directly: as explained above, the con- 
tours of constant li/Q become too closely spaced as they 
approach the neoclassical line in Fig. ITlT a). However, 
the temperature gradient of the final state, labelled B on 
both Fig. [6]and Fig. [TlTa.b). can be calculated indirectly 
by using the value of 7^; from Fig. [SJa) and rearranging 
(PSj) and ([57)1 to give k as a function of 7^, Pr^, W/WgB 
and Q/QgB- 



Vlt {V2qoRo/ro) jEQ/QgB 

n/Hgij + {iRlqo/vthiropDlE iXn 



(50) 



Taking — 1.17, Pr^ — 0.55 yields a temperature gra- 
dient at point B of K = 9.7. The temperature gradient at 
point A (the other side of the jump) is 7.4. Alternatively, 
using our model, the values of the temperature gradient 
can be read from Fig. [TlT b') which results in an identical 
jump oi K = 7.4 to K = 9.7. If we compare with the 
case of zero momentum input and zero flow shear (point 
C on Fig. [Ufa)), we find that flow shear has enabled 
a total increase in the temperature gradient of a factor 
9.7/4.5 — 2.2. If this were reproduced across the whole 
device, it could increase the core temperature by a fac- 
tor of around 10, but this would require a more detailed 
ID model (see^), rather than the OD model presented 
here. Note that while the jump in temperature gradient 
caused by the bifurcation (the transition from A to B) 
is a striking feature, a comparable contribution to the 
increase of the temperature gradient arises from the in- 
cremental suppression of the turbulent transport by the 
velocity shear (the difference between A and C). 



C. Neoclassical Heat Flux, Turbulent Momentum Flux 

It was noted earlier that the reduced transport state 
produced by the bifurcation is still turbulent, with a mo- 
mentum flux far greater than its neoclassical value (Fig. 
|6]). In Fig. [TlT b). it can also be seen that the reduced 
transport state docs not lie exactly on the neoclassical 
line; it docs however lie very close to it, which implies 
that the turbulent heat flux is much smaller than neo- 
classical. Thus, the bifurcation takes the system into a 
state where the heat is mostly transported neoclassically, 
hut the angular momentum is mostly transported by tur- 
bulence. The small ratio of Qt to Q„ reflects the fact 



Starting from point A, if Q were to be reduced at con- 
stant n/Q, the system would again be forced to jump 
to point B. In effect, what would happen is that the 
heat flux would become principally neoclassical; the mo- 
mentum flux would drop because Pr„ ^ Prt, and a bi- 
furcation would ensue in the manner described in Sec- 
tion IV El Thus, a decrease in the input heat could 
lead to a higher temperature gradient. We note, how- 
ever, that Q is normalised to the gyro-Bohm value: 
QgB = n-iTiVthipf/^V^RQ. Thus, decreasing Q/Qgs 
could correspond to decreasing the deposition of heat, 
but it could also correspond to an increase in tempera- 
ture or density^. 



E. Optimum Q 

Fig. [TlT b) shows that for every value of II/Q, there is 
an optimum Q that gives a maximum temperature gra- 
dient (the dotted line at k ~ 11 in Fig. [IHb)). The 
maximum temperature gradient increases with 11/ Q, but 
only slowly. The optimum value is the value of Q such 
that: 



dn 
dQ 



(51) 



n/Q 



and is plotted as a function of II/Q in Fig. [T^l 



F. Transition Region 

Finally, we will show that there is in fact only a lim- 
ited range of both Q and H/Q where bifurcations can 
happen in the way described above. To illustrate, we 
observe that no transition can occur along the contour 
n/Q = 0.12V2Ro/vth^ m Fig. [Ufb): if U/Q is kept con- 
stant at this value, decreasing Q / QgB from greater than 
its optimum value of 4.0 increases the temperature gra- 
dient smoothly up to its maximum value. The existence 
of a bounded region where such transitions can happen is 
studied in detail by the authors of^^, who calculate, using 
a simple transport model, the region in which transitions 
occur, and derive a criterion necessary for their existence. 
We apply the analysis using the model for the tur- 
bulent fluxes given in (|^5)) and (HHj), to determine the 
range of values of Q and H/Q for which bifurcations of 
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the form we have described can occur. In essence, bifur- 
cations can only occur when there exist multiple values 
of 7^; and k that give rise to the same values of Q and 
n/Q. From Fig. [TlTb). it is clear that this is only pos- 
sible for values of II/Q where there is a minimum in the 
curve of constant II/Q, and for values of Q which lie be- 
tween this minimum and cither a maximum in the curve 
or the point where the curve intercepts the neoclassical 
line. Thus the region in which transitions are possible is 
bounded by the curve 




10 100 1000 10000 100000 le+06 le+07 



dQ 
dn 



0, 



(52) 



n/Q 



and by the neoclassical line. This region is plotted in Fig. 

In order to give a clearer meaning to this diagram, we 
use equations pil - I33p and typical properties of plasma 
devices taken from the ITER Profile Database^ (and 
listed in Table |I| to replot the region in which transi- 
tions can happen in terms of the input beam power and 
beam particle energy, Pnbi and i^NBi • To generate these 
plots, we also calculate the collision frequency va self- 
consistently using p2)) . to take account of the varying 
strength of the neoclassical transport in each device. The 
transition regions for each device arc displayed in Fig. I13[ 
along with a point indicating the typical values of Pnbi 
and i?NBi in each device. Fig. [13] shows that in which 
transitions can happen lie within an order of magnitude 
of the typical values. It should be stressed that with the 
simplified magnetic geometry (Section |Tl| and model of 
neoclassical transport (Section IV Dp used in this study, 
and with the many assumptions about the way the en- 
ergy and momentum are deposited in the plasma (Sec- 
tion |VB]), closer agreement could not be expected. In 
particular, the assumptions of Section IV Bl are likely to 
overestimate the applied torque and hence overestimate 
the beam energy needed for a transition. 




0.02 0.04 0.06 0.08 0.1 



0.12 0.14 



FIG. 13. The regions (shaded) in wliich transitions can hap- 
pen, as a function of the beam power and the beam particle 
energy, for three different devices. The dashed lines indicate, 
for each device, the value of Pnbi for a given Bnbi which 
would lead to the optimum temperature gradient, as described 
in Section [VII El The points indicate typical values of Pnbi 
and Pnbi for each device. The projected Pnbi and Pnbi for 
MAST Upgrade were taken from Ref.— . 



TABLE I. Typical plasma properties from the ITER Profile 
Database^. The symbol a denotes the minor radius of the 
device. The temperature was calculated as the mean of the 
core (TIO) and edge (TI95) temperatures. 

Ti {eV) Ui {m-^) R (m) Bt (T) a (m) 

DIII-D 2.7e+03 9.0e-hl9 1.7e+00 1.7e-F00 6.1e-01 

JET 2.6e+03 6.3e-hl9 2.9e-h00 2Ae+00 9.3e-01 

MAST 0.6e+03 3.9e-|-19 8.0e-01 5.3e-01 5.6e-01 



VIM. CONCLUSIONS 

We have determined the way the turbulent transport of 
heat and toroidal angular momentum is affected by the 
equilibrium gradients of temperature and velocity over 
a wide range of those two parameterSjin the case of a 
plasma with a sonic background flo-wrH. We have ex- 
tended the range of flow shears up to 20 times the linear 
ITG instability growth rate calculated at zero flow shear. 
We have considered the zcro-magnetic-shear regime. We 
have used a low value of the safety factor, go = 1-4, which 
reduces the strength of the destabilising parallel velocity 
gradient drive relative to the stabilising perpendicular 
velocity shear. In so doing we have established the fol- 
lowing key results: 

1. The plasma is linearly stable to microinstabilities 



FIG. 12. The region in which bifurcations can occur, calcu- 
lated using the analysis o^ and the parameterisation of the 
fluxes described in Section IVTl The cross indicates the loca- 
tion of the bifurcation described in Section IV El The dashed 
line represents the optimum value of Q for a given value of 
n/Q, Eq. JSU. 



We have neglected the effect of the Coriolis force due to such 
a flow in our study. However, previous work^ indicates that its 
effect is to reduce the effective turbulent momentum diffusivity 
through a momentum pinch. This will make higher-flow-shear 
regimes easier to access but will not affect the bifurcation mech- 
anism, since the neoclassical transport of momentum remains 
much lower. 
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for all non-zero values of flow shear but subcriti- 
cal turbulence can be sustained nonlinearly across 
much of this stable region. 

2. For velocity shear 7^; ^ 1 (we remind our readers 
that is normalised to Vthi/Ro), increasing 7^ 
reduces the levels of turbulence driven by the ion 
temperature gradient. For 7^ > 1, the destabilis- 
ing effect of the parallel velocity gradient becomes 
dominant and increasing the flow shear increases 
the levels of turbulence and turbulent transport. 

3. For low flow shears, -fE < 0.8, and low heat flux, 
Qt/QgB ^ 2.5, flow shear reduces the thermal 
transport stiffness. For low flow shears and higher 
Qt, flow shear increases the thermal transport stiff- 
ness. At high flow shears, 7b ^ 1, the thermal 
stiffness is reduced as the PVG drive becomes dom- 
inant. 

4. For 7b 0.1 the turbulent Prandtl number stays 
within the range 0.5-0.8 across a wide range of 
temperature gradients. Thus, the turbulence trans- 
ports heat and momentum with comparable vigour. 

5. For Rq/Lt ^ 11, there is a large region around 
7b ^ 1 where the turbulence is completely 
quenched. 

6. Owing to the existence of this region, and the fact 
that the neoclassical Prandtl number is much lower 
than the turbulent Prandtl number, it is possible 
to trigger a bifurcation to a regime of high velocity 
and temperature gradients by either (i) increasing 
the momentum input at constant heat input or (ii) 
lowering Q/Qgs at a constant ratio of heat to mo- 
mentum input. 

7. In the high gradient (reduced transport) state that 
is reached via this bifurcation, the heat is princi- 
pally transported neoclassically, whereas the mo- 
mentum is principally transported by turbulence. 

8. For any given input of toroidal angular momen- 
tum, there is an optimum input of of heat which 
maximises the temperature gradient. Increasing 
the input of heat above this optimum can reduce 
the temperature gradient by increasing the levels 
of turbulence. 
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